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Abstract 



We study the KPZ equation (in D = 2, 3 and 4 spatial dimensions) by using a RSOS 
ly-) ■ discretization of the surface. We measure the critical exponents very precisely, and we show 

that the rational guess is not appropriate, and that 4L> is not the upper critical dimension. We 
are also able to determine very precisely the exponent of the sub-leading scaling corrections, 
that turns out to be close to 1 in all cases. We introduce and use a multi- surface coding 
|Q ! technique, that allow a gain of order 30 over usual numerical simulations. 

C3 



1 Introduction 



The KPZ equation |I|], in its apparent simplicity, involves many issues that need clarification. The 
O continuum equation (that describes the local growth of an interface profile) is 

^l = u V 2 h + ^(Vh) 2 + i 1 (r,t) . (1) 

The real behavior of this equation is not well understood. The absence of a complete mean field 
theory does not help, and the fact that we have do understand the effects of a strong coupling fixed 
point makes very difficult the development of a perturbative renormalization approach. 

A practical approach for numerical studied of the problem is to consider lattice Restricted Solid 
On Solid (RSOS) surfaces (see for example for a review): one considers a discretized height field 
on a D dimensional lattice, and imposes the constrain that the absolute value of the distance among 
two neighboring surface elements can only take the values zero and one. 

Old extensive numerical simulations 0, |J do not clarify the situation completely. After a number 
of interesting Field Theoretical results [f|, the introduction of new Renormalization Group based 
techniques is a potentially promising direction of development [H. 
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The main theoretical points which still deserve proper explanation are two: it is still not clear 
whether equation ([[]) has a finite upper critical dimension (.D>) or not, and the exact quantification 
of the related critical exponents. Regarding arguments in favor of D > = 4 we address the reader to 
0, while arguments supporting D > = oo can be found in || |8|]. 

In this light we introduce here a new numerical technique and present some precise numerical 
simulations that allow us to estimate critical exponents. Thanks to the new, precise technique we 
succeed to clarify some questions: we show for example that a conjecture of || is not founded, and 
we give quantitative estimates about the behavior of sub-leading corrections. 

2 The Multi-Surface Coding 

The precise results on which this paper is based are due to a new technique we introduce to simulate 
RSOS surfaces, where adjacent elements of the surface can only be at a distance of 1, or — 1 
lattice spacings. The technique is a generalization of the so called multi-spin coding technique (well 
discussed by Rieger |T(J in the context of the study of disordered systems), where, by using the fact 
that the ±1 spins can be represented by Boolean logical variables, one stores 64 copies of the system 
in one single computer word (we will assume in the following we are using a 64 bit computer) . The 
main idea is that we can simulate, basically at the cost of one single usual simulation, 64 copies 
of the system, by rewriting the basic operations (like summing spins for computing the effective 
force) as boolean operations, and by exploiting the fact that when, for example, the computer is 
calculating a logical bit AND it is indeed doing that 64 times at once. The gain of such an approach 
over an usual single spin simulation is of order 30: one gains a factor 64 for the number of systems 
that updates at once, and looses a factor of order 2 in computational complexity (adding the spin 
in the Boolean form takes some more time [|10|j ). 

We generalize here the Ising case to the case of a field of differences, that can have three values, 
since a given element of the surface can be at the same height of a given neighbor, or one step 
behind it, or one step ahead of it. The method is new as applied to such a system, where one does 
not have to compute the value of an energy, but uses the boolean operations to determine if the 
element can be moved without violating the geometrical constraint. 

Let us thinks for simplicity about 2D, where we have a 2D support where a height field (the 
surface height) hi = h X)V is defined. One surface element hi has 2D = 4 first neighbors: in the RSOS 
model the difference A ijjU = hi — hi±^, where /i = x,y, can only take the three values —1, and +1. 
We store the DL 2 values of the Aj i(U (and, as we will discuss better, each A i)M needs two bit to be 
stored. This is a redundant way to store the information, since a factor two in memory could be 
easily saved, but it is very convenient from the point of view of computer time: as usual one trades 
memory for time, and avoids a cumbersome reconstruction by storing more information). 

If a given element is behind its neighbor the difference is —1, if it has the same height than 
the neighbor the difference is 0, while if it is ahead of the neighbor the difference is +1. Since we 
have three allowed values we can represent each of these difference with two bits (that we will call 
respectively H, high bit, and L, low bit): we can for example code with 00 the situation where the 
given element is behind its neighbor, with 01 the situation where they have the same height, and 
with 10 the situation where it is ahead (the value 11 is forbidden). 

If the element is ahead of even a single one of the 2D neighbors the move is forbidden (it would 
violate the RSOS constrain). It is clear this is easily implemented in our coding: one just needs to 
do a logical OR of the 2D H bits related to the site i, and if it is 1 at least one of the neighbors 
is behind and the move is not allowed. Let us see better why. We start comparing our element i 
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to the first of its 2D = 4 neighbors (let us say the one in the positive direction +1): if i is ahead 
the neighbor the relevant H bit (that we have stored in A ij+1 ) is 1, and we already know we cannot 
move. If i is on the contrary at the same height of the neighbor or behind it we have that, as far 
as this neighbor is concerned, the element % can advance of one unit. Now we look at the second 
neighbor, where the same reasoning holds: if we look at the logical OR of the two relevant high 
bit we will find that we cannot move if this quantity is one. Looking at all the 2D neighbor we 
see that the move is forbidden if the logical OR of the 2D H bit is 1: clearly this operation, as all 
the other present in the core of the code, updates with a single computer cycle the 64 copies of the 
system. After doing that, if the element i cannot be moved there is nothing that has be done: if it 
gets moved we have to update the A ijjU related to the site % and to the neighbors, to describe the 
new situation. 

The codes simulating systems in a different number of spatial dimensions (in our case from 
D = 2 to D = 4) are simple generalization one of the other (when adding a dimension one has to 
add checks on the new neighbors and their update). 

As usual in these kind of simulations the random noise is implemented with a random choice 
of the site to be updated. We also implement a fifty percent probability of really updating a 
surface element that according to the RSOS constrain could be updated. This is very important, 
since starting from random independent surfaces is not enough: because of our parallel scheme 
the sites of our 64 copies have to be updated in the same order, and such an updating algorithm 
with updating probability one is attractive, and the 64 configurations asymptotically at large times 
become equal [|ll|]. The probability of not accepting an allowing change, that depends on each of 
the 64 configurations, solve this problem. 

We are aware of a parallel algorithm to update surfaces [Q: it is very different in spirit from 
our algorithm (since it parallelizes on different sites of the lattice). We have not compared in detail 
the performances of the two algorithms, but we believe that on one side the algorithm of [0 is 
more general, and not limited to RSOS models, but on the other side our algorithms is more regular 
(there are no exceptional loops), and is probably better performing for the model we study. 

3 The Numerical Simulation 

We have based this work on numerical simulations of lattices of volume V = L D . The spatial 
dimensionality D is the spatial support, where a 1 dimensional surface take values. We study the 
D = 2, D = 3 and D = 4 cases. At a given time t of the dynamical evolution the position of the 
surface can be expressed by the values hi(t) (that we reconstruct by the differences we store in our 
code, see the former section), i represents, in lexicografic order, a D-ple labeling the spatial sites. 

We consider a dynamics which generates a Restricted Solid on Solid (RSOS) growth: distances 
of first neighboring elements of the surface cannot be larger than one. At each trial step we move 
elements of the surface that are not constrained not to do so because of the RSOS restriction with 
probability |: we have to use a probability different from one in order to keep independent the 
different surfaces we simulate in the same computer word. 

We consider a large number of different lattice sizes. In D = 2 we take L going from 5 to 641, 
in D = 3 we consider L going from 5 to 103 while in D = 4 L goes from 5 to 28. All our data have 
been averaged over 64 different dynamical runs. 

Because of the way we use to determine critical exponents, by trying to measure precisely the 
asymptotic time behavior, we use very long runs, and we always try to check that we have reached 
the asymptotic plateau in a clear way (see the discussion and the figures in the next section). Let 
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Table 1: Number of full lattice sweeps for each run on different lattice sizes and number of 
dimensions. 



t be the time labeling sweeps of our simulation (we define a sweep as the trial update of V random 
sites). In any of our simulations we run T updating sweeps: we give in table p] the number of full 
lattice sweeps for each run on different lattice sizes and number of dimensions. We measure the 
observable every 1000 lattice sweeps: when analyzing the large time asymptotic behavior of the 
system we discard the first half sweeps (we call T the time of our first measurement): this is a very 
conservative attitude, but we prefer to be safe on not having any systematic bias by paying a price 
of making maybe the statistical error ten percent larger of the best we could do. 
We define the time dependent observables 



= 77 EN*)] 

i=l 



(2) 



and 



h(t)-h(t) 



(3) 



that we compute for k = 2, 3 and 4 and for different D and L values (when needed we will label w 
with the upper script (L), to make clear to which lattice size we are referring). We define the large 
time asymptotic limit of Wk(t) as 



1 - 1 + 1 t=To 



(4) 



for simulation on a lattice of linear size L. We always check (and this is one of the crucial points of 
this note, that we will discuss in better detail in the next section) that T and T are large enough 
to make our result unbiased in the precision of our statistical error. 



4 Analysis 

We will discuss here the analysis of our numerical data. We want to determine the critical exponents 
of the asymptotic behavior of the Wk we have defined in the former section. For example we have 
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that the asymptotic infinite time value of w 2 has a leading scaling behavior 



wi L) = L 2x , (5) 

while at intermediate times (large enough for being in the scaling region but small enough not to 
feel the finite size of the lattice) 

w 2 {t) = t~ . (6) 

The first behavior is obtained by taking large times on different lattice sizes, and by studying the 
time asymptotic value as a function of L. The second behavior is studied by simulating large 
lattices, and analyzing the behavior of the systems for times larger than one, but very smaller than 
the thermalization time (at the given value of the lattice size). An exact (Galilean) invariance of 
the KPZ equation implies that z + x — 2- Deciding if it is better to measure accurately z or \ is a 
practical matter. 

As opposed to the choice, for example, of reference ||, here we have mainly based our analysis 
on fitting the behavior of the large time asymptotic value w[ L \ and we have only used fits to the 
intermediate time behavior to substantiate our results. We believe in fact that it is very difficult to 
determine a precise quantitative estimate of the exponent of the time scaling. The problem with 
the time dependent behavior of Wk(i) is that, in order to get an unbiased value, one needs a double 
cutoff, both at small and at large time. At small times the behavior of Wk(t) is not a pure power, 
and one has to discard small lattice, and/or to use corrections to scaling, in order to remove this 
effect. At large times one starts to feel the finiteness of the lattice system, and a new crossover 
(toward the asymptotic, constant time behavior) intervenes. In other words a careful analysis of 
the time exponent needs a double sliding window, moving both at small and at large time. We also 
find that the crossover effects at large time are very important: even on large lattices one can soon 
see systematic effects on the exponent estimate due to the finiteness of the lattice. 

On the contrary the time asymptotic behavior only needs one cutoff, that excludes small times, 
where the asymptotic value has not yet been reached. This is easy, and we do it by using a 
logarithmic division of our data. We are in this way able to check with very high precision that 
we are computing an unbiased (effective, L-dependent) exponent. Again, we will also show results 
obtained from a direct fit of z, to show they are consistent with the \ values we determine. 



Our main analysis is done by fitting at the same time the three sets of data [T3| : 



w 2 ~ A 2 L 2x (1 + B 2 L~ 

w 3 ~ -A 3 L 3 * (l + B 3 L-") , (7) 
~ AaL 4x (i + BaL 



We always compare this fit to the best fit of w 2 alone (typically by only including the large volumes 
and by ignoring scaling corrections) and check that things are coherent. We have also check inde- 

3 

pendently that, for example, w 3 really scales as w 2 : this is clearly true for our data. Without the 
use of all our set of data that we have defined in ([/]) we would not have been able to determine u 
with a reasonable statistical precision. 

With this definition the exponent % ls t ne same of \ of Q and of a of |6j (our definition of 
dimensionality of the system excludes the time dimension, and is always only the dimension of the 
space). 



5 



as 




-0.265 -0.26 -0.255 -0.25 -0.245 -0.24 



Figure 1: _R 4 versus R 3 , D = 2. 



The error analysis is done by using a jack-knife approach [13]: we divide our statistical sample 
in 10 parts all including all of data but one tenth (each part excludes a different tenth of the data), 
we fit ten time the behavior of, for example, (|7[), and compute the error on Ak, -Bfc, X an d oj by 
the fluctuations of the ten results (multiplying the error times a factor 10, due to the fact that the 
individual parts we have formed are correlated ||14j| ). In short, we present a reliable estimate of the 
statistical errors over the quantities we determine. 

The first point we have to discuss is about the exponents of finite size corrections that appear 
in equation (|7|). Do the corrections to even and odd momenta really scale with the same exponent? 
This has been questioned in [15], and we provide here accurate evidence that U2 = ^3 = uo^ = u> in 
our model. 

Let us start from arguing that we are not in a situation in which the scaling exponent of the odd 
moments, u Q dd is smaller than the scaling exponent of the even moments, uj even (this is the opposite 
scenario of the one proposed in |T5|). 



Following [15 1 we define 



w 3 



Ra 



U>4 



The exponent x disappears from these ratios, and, on general ground, if u dd < u e 
asymptotically for large L, ignoring sub-leading corrections, 



one has that 



R 3 ~ c 3 + d 3 r° d 



R A ~ c 4 + di 



(9) 



i.e. 



R 4 ~ c + d (R3 + c ) uleven/u} ° dd , (10) 

and i?3 is a linear function of R^ if and only if 0J o dd > ^even (i n which case the two ratios asymptoti- 
cally scale with the same exponent). We plot in figure [I] R4 versus R3, and notice that the linearity 
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Figure 2: L dlog ^ (L)) versus \ in 2D, for k = 2,3, 4. 



is impressive. In these data there is no sign of a discrepancy among scaling exponents of odd and 
even momenta, and they surely exclude that u dd < u even . The last point on the left in the figure is 



our asymptotic extrapolation, with attached the (small) statistical error (the best estimate of (15 



is i?3 = —0.27 ± 0.01 and i?4 = 3.15 ± 0.02, well compatible with our data but with a very larger 
error). 

We use now figure || to also exclude the case u dd > u even , establishing in this way that for our 
model u dd = ^even- m figure |^ we plot the effective scaling exponent obtained by using separately 
the data for w 2 (L), w 3 (L) and u> 4 (L). The three quantities do all depend linearly, with very good 
approximation, on j, showing that we are having the same exponent of the scaling corrections. 
Again, the impressive linearity of the data implies that we are measuring precisely a single exponent 
of finite size corrections, and also that we are not mislead by finite size corrections. 

Let us now discuss the determination of the exponent \. 

Let us start by discussing the D = 2 case. A simple analysis of w 2 without corrections to 
scaling shows that a fit to lattices of linear size from L = 19 give a good value of the chi squared 
and x — 0.393. A systematic analysis to the form (0) by including lattice with L > 11 gives our 
final best value of 



XD=2 = 0.393 ± 0.003 , u> D =2 = 1.1 ±0.3 . (11) 

We plot the rescaled w 2 versus the rescaled time in figure ||]: it is clear that the asymptotic plateau 
is exposed with good accuracy, and that the scaling is very good. Also the best fit to the form (|7|) 
is very good: we plot the numerical data for w 2 , w 3 and u> 4 versus L and the best fit in [|. 

We can compare to the rational guess of that would give here \R — 0.4. Indeed Kim 
and Kosterlitz in || conjecture that x(d) = ^3 (that seemed to fit reasonably the numerical 
results available at the time). Here Lassig |16|], by using an operator product expansion, also find 
x{d = 2) = 0.4. Our result is at three standard deviations from the rational guess, that is a safe 
distance. Still, since we are dealing with a very complex situation, with many corrections that can 
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Figure 3: Scaling plot of the rescaled w 2 versus the rescaled time. 2D. 
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Figure 5: As in figure |, but 3D. 

possibly affect the result (sub-sub-leading corrections, short time, small volume,...), we perform a 
further check to determine if x/? = 0.4 is a plausible result. We fix x — 0.4, and fit our data with 
now seven and not eight free parameters (A 2 , A3, A4, B 2 , B 3 , B4, and oS). Now we get a very small 
value of u ~ .28, and a chi squared that increases of a factor 10 from our previous best fit (where 
it was of order one per degree of freedom). 

In order to show that the fact that we have been able to exclude that the exponent takes the 
value 0.4 is not due to the hypothesis that the exponent of the sub- leading correction is the same 
for all momenta we can look again at figure |^. Already from these data it is clear that the value 0.4 
is excluded. The more sophisticated analysis which we have presented before is crucial in obtaining 
a controlled extrapolation to L = 00, keeping the statistical errors under control. We believe that 
this is a very strong evidence against the validity of the rational guess. We will see that for higher 
D values we get an even clearer discrepancy. 

In 3D the same analysis of the three cumulant allows to establish that 

XD=3 = 0.3135 ±0.0015 , u D=3 = 0.98 ± 0.08 . (12) 

We are including here sizes from L = 11 up to L = 103. We plot the rescaled w 2 versus the rescaled 
time in figure ^| and the numerical data for u> 2 , ^3 and 104 versus L and the best fit in ^. Here the 
rational guess would give Xr — § — 0.333, and we are sitting at more than ten standard deviations 
from it (Lassig gives here | that is far from our estimate). The same check we have done in D = 2 
leads to a strong evidence: when fixing x = | we find again a very small value u ~ 0.2, and the 
chi squared increases of a factor large than 20. Here the evidence against the validity of a rational 
guess is even stronger than in D = 2. 

In 4D we use data for L starting from 10. Again the same three cumulant analysis gives us 

Xd=a = 0.255 ± 0.003 , u) D=i = 0.98 ± 0.09 . (13) 
We plot the rescaled w 2 versus the rescaled time in figure ^ and the numerical data for w 2 , w 3 and 
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Table 2: Best fit values of the parameters entering equation (0). 



W4 versus L and the best fit in |8|. 

Here xr — f — 0.286, and, again, we are sitting a ten standard deviations away from the rational 
guess. 

We summarize our findings in table where we give all the best fit values of the parameters 
entering equation (|7|). 

Two more observations are interesting. In first we find uj ~ 1 for all the D values we investigate. 
In second the pre-factor of the scaling corrections increases with D: we find B 2 (D — 2) ~ 0.08, 
B 2 (D = 3) ~ 0.25, B 2 (D = 4) ~ 0.37. 



Our recent are not incompatible with the recent ones of | 1yB , but our small error bars allow us 
to reach precise conclusions. Also the comparison with the exponents found in 0] is fair: we have 
stressed the length of our runs, in order to be able to give a clean estimate of the time asymptotic 
behavior, so that our result is hopefully unbiased. 



5 Conclusions 

The numerical technique we have introduced works well, and has allowed us to run very precise 
numerical simulations with a limited amount of computer time (a few months of Pentium II pro- 
cessor) . 

We have been able to determine critical exponents of the KPZ universality class with high 
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accuracy. We have falsified the guess that the exponents are simple rational numbers. It is also 
unambiguous from our data that the upper critical dimension is larger than 4 (as opposed to the 
claims of ]7j). 

Thanks to our precise measurements (and fitting together the first 3 non trivial moments of h) 
we have also been able to determine the exponent u of the first non-leading scaling corrections. It 
is interesting to notice that the estimated value of u is always very close to 1, independent from the 
dimensionality of the system. 

The next interesting step, following the approach of ||, would be to try and implement a system- 
atic Monte Carlo Renormalization Group: the Mult i- Surface Coding technique we have discussed 
could be a very important ingredient of such a development. 
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